🧮 Size-selection of diatoms in fluctuating environments - Quantitative Approaches to Organismal Biology

This notebook implements (with some modifications) the diatom population dynamics model presented in Litchman et al..

This model uses allometric trends in growth rates and storage capacities to calculate the population dynamics of different-sized diatom variants. The variants compete for a nutrient (nitrogen in the form of nitrate) within a mixed layer. Nutrient is depleted in the mixed layer due to uptake by diatoms.

Mixing events (e.g. storms) occur periodically, that replace a fraction of the water in the mixed layer with deep water. This has two effects: It transports some of the diatom cells into deep water, where they ultimately die; and, it injects nutrient (contained in deep water) into the mixed layer. Therefore, a mixing event both lowers the diatom population and provides nutrients to enhance its growth.

Litchman et al.‘s model is a so-called quota model. “Quota” refers to cells’ storage of nutrients, increasing as cells take up nutrient from the water and decreasing as that nutrient is consumed during growth. The mode tracks the changes over time in:

The number of cells of each diatom variant
The quota (amount of stored nutrients) in cells of each diatom variant
The amount of nutrient in the water

The key question asked by Litchman et al. was,

Do variations in the frequency and intensity of mixing events, combined with allometries of nutrient storage and use, explain the size distributions of diatoms in different habitats?

The rationale and key parameters are explained in this discussion of diatom allometry and this discussion of Litchman et al.'s data analysis. Explanation of the formulation of this model is in this context page (and in the original paper).

How to use this notebook¶

To use the model, start by selecting Run All Cells under the Run menu at the top left of this page.

The parameters used in simulating diatom population dynamics are presented in three sets of textboxes and execution buttons, which you will use in sequence to:

Set up the number and size range of diatom variants, and the characteristics of the mixed layer they inhabit.
Seed the mixed layer with one or more of the diatom variants.
Set the interval between mixing events and how many days to simulate the population.
Run the simulation and look at the output.
Optionally, reseed the population with new cells, and restart the simulation to see how the new set of variants interact.

From the initial conditions you specify, the model runs through many mixing intervals and the diatom populations approach a repeated pattern.

Graphical output¶

The model produces two sets of graphical output:

Plots of the population (expressed as biovolume) and quota of each variant, and the nutrient concentration, through the last simulated mixing interval.
These “period” plots are good indicators of how effectively each variant stores nutrients when it is available, and uses that storage to fuel growth.
Plots of the averages over each mixing interval of the population and quota of each variant, and the nutrient concentration, over the entire simulation.
The averaged plots are good indicators of the outcomes of competition between variants, because they show positive or negative population trends over long timespans.

Biovolume is the number of cells times the volume of each cell, for each of the size variants. If cell densities are approximately constant, then the biovolume of each variant is approximately proportional to the total biomass of that variant.

These plots are automatically saved at the end of each simulation as images, so you can copy and paste them into a document or presentation to record your results. The png format images are period.png and averagred.png; the svg format images are period.svg and averaged.svg.

Tabular output¶

The model also produces two tables of period-averaged results for the final mixing period:

A table of standing stock, averaged through the last mixing period.
This table shows, for each size class:
- class number ( $i$ )
- size metric ( $s_i$ )
- size ( $10^{s_i}$ )
- number cells per liter ( $N_i$ )
- relative quota ( $\frac{Q_i}{Q_{i,min}}$ )
- biovolume ( $BV_i = N_i \times 10^{s_i}$
- nutrient concentration, $R$
Note that $R$ reflects the entire mixed layer, and does not vary by diatom size – it is included in each row only to make a complete table for import into a spreadsheet.
A table of diatom population losses (in units of biovolume per day), averaged through the last mixing period.
This table shows, for each size class:
- class number ( $i$ )
- size metric ( $s_i$ )
- size ( $10^{s_i}$ )
- losses due “background” mortality, $M$
- losses due to sinking out of the mixed layer, $S$
- losses transport to deep water by mixing events ( $L$ )

# Load modules and set graphics environment
%matplotlib widget
import numpy as np
from matplotlib import pyplot as plt
#plt.ion();
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets
from math import *
from Litchman2009 import diatoms

# Instantiate a Params object and a GUI to modify parameters. 
#params=Params(auto_plot=True)
h = '27px'
n_sizes=widgets.IntText(value=8,width=10,description = r"$n_{sizes}$")
L_n_sizes = widgets.Label(value='Number size classes',layout=widgets.Layout(display="flex", 
                                        justify_content="flex-end",width="100%",height=h))
min_size=widgets.FloatText(value=2.5,
                               description = r"$S_{min} (log_{10})$")
L_min_size = widgets.Label(value='Min. cell size',layout=widgets.Layout(display="flex", 
                                        justify_content="flex-end",width="100%",height=h))
max_size=widgets.FloatText(value=2.5+1.25*8,
                               description = r"$S_{max} (log_{10})$")
L_max_size = widgets.Label(value='Max cell size',layout=widgets.Layout(display="flex", 
                                        justify_content="flex-end",width="100%",height=h))
qinit=widgets.FloatText(value=0.1,description = r"$q_{init}$")
L_qinit = widgets.Label(value='Initial fractional quota',layout=widgets.Layout(display="flex", 
                                        justify_content="flex-start",width="100%",height=h))
m=widgets.FloatText(value=0.08,description = r"$m$")
L_m = widgets.Label(value='Mortality rate ($day^{-1}$)',layout=widgets.Layout(display="flex", 
                                        justify_content="flex-start",width="100%",height=h))
z_m=widgets.FloatText(value=25,description = r"$z_m$")
L_z_m = widgets.Label(value='Mixed layer depth ($m$)',layout=widgets.Layout(display="flex", 
                                        justify_content="flex-start",width="100%",height=h))
a=widgets.FloatText(value=0.3,description = r"$a$")
L_a = widgets.Label(value='Fractional mixing',layout=widgets.Layout(display="flex", 
                                        justify_content="flex-start",width="100%",height=h))
Rdeep=widgets.FloatText(value=40.,description = r"$R_{deep}$")
L_Rdeep = widgets.Label(value='Deep nitrogen',layout=widgets.Layout(display="flex", 
                                        justify_content="flex-start",width="100%",height=h))
ui1 = widgets.VBox([n_sizes,min_size,max_size,m])
ui2 = widgets.VBox([a,Rdeep,z_m])
uiL1 = widgets.VBox([L_n_sizes,L_min_size,L_max_size,L_m])
uiL2 = widgets.VBox([L_a,L_Rdeep,L_z_m])

Setting up diatom variants and mixed layer characteristics¶

This part of the notebook lets you set up the “scenario” of diatom variants and mixed layer characteristics.

The parameters are:

$N_{sizes}$ , the number of diatom variants
$S_{min}$ and $S_{max}$ , the minumum and maximum sizes of the variants (variants will be distributed across this range)
$m$ , the mortality rate due to predation, pathogens, etc., affecting all variants equally
$z_m$ , the depth of the mixed layer; a shallower layer means more loss due to sinking, which disproportionately affects larger cells.
$a$ , the fraction (between 0 and 1) of the mixed layer that is replaced in each mixing events. Larger $a$ means more cells are lost, but more nutrients are gained.
$R_{deep}$ , the nutrient concentration in deep water.

It is a good strategy to start with these default values, to see their results, then modify one of more to see the effects of those changes.

When you have the parameters you want, click the Set up population button below to set up the mixed layer.

# Set up textboxes for population characteristics
isizes=widgets.Text(value='all',description = r"$i$")
L_isizes = widgets.Label(value='Size classes ($S_i$) to seed',layout=widgets.Layout(display="flex", 
                                        justify_content="flex-start",width="100%",height=h))
Nseed=widgets.FloatText(value=10,description = r"$N_i$")
L_Nseed = widgets.Label(value='Number cells to seed',layout=widgets.Layout(display="flex", 
                                        justify_content="flex-start",width="100%",height=h))
t_mix=widgets.FloatText(value=8.,description = r"$t_{mix}$")
L_t_mix = widgets.Label(value='Mixing interval ($days$)',layout=widgets.Layout(display="flex", 
                                        justify_content="flex-start",width="100%",height=h))
n_pers=widgets.IntText(value=32,description = r"$n_{pers}$")
L_n_pers = widgets.Label(value='Number of intervals',layout=widgets.Layout(display="flex", 
                                        justify_content="flex-start",width="100%",height=h))
n_record=widgets.IntText(value=12,description = r"$n_{record}$")
L_n_record = widgets.Label(value='Data per interval',layout=widgets.Layout(display="flex", 
                                        justify_content="flex-start",width="100%",height=h))
ui3 = widgets.VBox([qinit,isizes,Nseed])
ui4 = widgets.VBox([t_mix,n_pers,n_record])
uiL3 = widgets.VBox([L_qinit,L_isizes,L_Nseed])
uiL4 = widgets.VBox([L_t_mix,L_n_pers,L_n_record])

uis12 = widgets.HBox([uiL1,ui1,ui2,uiL2])
uis3 = widgets.HBox([uiL3,ui3])
uis4 = widgets.HBox([ui4,uiL4])

global n_sizes_,min_size_,max_size_,m_,a_,Rdeep_,z_m_
def setup_pars(n_sizes,min_size,max_size,m,a,Rdeep,z_m):
    global n_sizes_,min_size_,max_size_,m_,a_,Rdeep_,z_m_
    # Put current parameters into global equivalents for later use
    n_sizes_ = n_sizes
    min_size_ = min_size
    max_size_ = max_size
    m_ = m
    a_ = a
    Rdeep_ = Rdeep
    z_m_ = z_m

out12 = widgets.interactive_output(setup_pars,{'n_sizes':n_sizes,'min_size':min_size,'max_size':max_size,
                                          'm':m,'a':a,'Rdeep':Rdeep,'z_m':z_m})
display(uis12,out12)

# Set up a new diatoms object, using current parameters
global D
buttonD = widgets.Button(
    value=False,
    description='Set up population',
    disabled=False,
    button_style='', # 'success', 'info', 'warning', 'danger' or ''
    tooltip='Click to set up a simulation model of a diatom population, using the current parameters.',
    icon='' # (FontAwesome names without the `fa-` prefix)
)
outputD = widgets.Output()
@outputD.capture()
def setup(b):
    global D, n_sizes_,min_size_,max_size_,m_,a_,Rdeep_,z_m_
    # Initialize diatoms object
    D = diatoms(z_m=z_m_,n_sizes=n_sizes_,size_range=[min_size_,max_size_],
                m=m_,a=a_,Rdeep=Rdeep_,static_v=False)
    D.setup()
    print('New diatom population initialized...')
buttonD.on_click(setup)
widgets.VBox([buttonD,outputD])

Seeding the mixed layer with diatom variants¶

This part of the notebook lets you “seed” diatoms by adding cells into the scenario, as a starting point for simulating population dynamics. When you set up the mixed layer, the populations of all variants were set to zero. Using these textboxes, you can introduce a starting population of diatoms, either all variants simultaneously or only variants you specify.

The parameters are:

$q_{init}$ , the initial quota in new diatom cells at the start of the simulation
$S_{i}$ , the size classes of the variants with which to seed the mixed layer.
- If this entry is “all”, then all variants will be seeded
- If this entry is one or more numbers, e.g. “1” or “2,6” then only those size classes will be seeded
$N_i$ , the number of cells to seed into each selected size class.

It is a good strategy to start with these default values, to see their results, then modify one of more to see the effects of those changes.

When you have the parameters you want, click the Seed cell population button below to set up the scenario.

# Generate textboxes for parameters to seed the population
global qinit_,isizes_,Nseed_
def seed_pars(qinit,isizes,Nseed):
    global qinit_,isizes_,Nseed_
    # Place parameters for seeding cells into global variables for later use.
    # This is executed on each parameter change.
    qinit_ = qinit
    isizes_ = isizes
    Nseed_ = Nseed

out3 = widgets.interactive_output(seed_pars,{'qinit':qinit,'isizes':isizes,'Nseed':Nseed})
display(uis3,out3)

buttonS = widgets.Button(
    value=False,
    description='Seed cell population',
    disabled=False,
    button_style='', # 'success', 'info', 'warning', 'danger' or ''
    tooltip='Click to add cells to the current population using the current parameters.',
    icon='' # (FontAwesome names without the `fa-` prefix)
)
outputS = widgets.Output()
@outputS.capture()
def seed(b):
    global D,qint_,isizes_,Nseed_
    # Seed selected size class(es), using the parameters set in seed_pars.
    # This is executed only when the button is clicked. 
    # isizes=='all' means all size classes
    jsizes = eval(isizes_)
    if isizes_ == 'all':
        D.seed(qinit=qinit_,N=Nseed_,sizes='all')
        print(f'Seeded all size classes with {Nseed_} cells per liter')
    elif isinstance(eval(isizes_),tuple):  # multiple integers entered
        if any(i >len(D.sizes)-1 for i in eval(isizes_)):
            print(f'Size class is outside range -- please choose a different size class.')
            return
        D.seed(qinit=qinit_,sizes=list(eval(isizes_)),N=Nseed_)
        print(f'Seeded selected size classes: {isizes_} with {Nseed_} cells per liter')
    elif isinstance(jsizes,int):   # a single integer entered
        if jsizes>len(D.sizes)-1:
            print(f'Size class is outside range -- please choose a different size class.')
            return
        #sizes_ = jsizes #[].append(jsizes)
        D.seed(qinit=qinit_,sizes=jsizes,N=Nseed_)
        print(f'Seeded selected size classes: {isizes_} with {Nseed_} cells per liter')
    else:   # none of the above -- ask for new input
        print(f'Invalid size selection -- please choose a different size class.')
        return
    print(f'Current cell populations are {D.Ns}')
buttonS.on_click(seed)
widgets.VBox([buttonS,outputS])

Setting mixing event intervals, and running the simulation¶

This part of the notebook lets you specify the interval between mixing events, to see what the effects different intervals have on competition between diatom variants.

Important: These parameters run the simulation using the existing diatom populations. If you want to start with a “blank slate”, go back and click the Set up population and Seed cell population buttons.

The parameters are:

$t_{mix}$ , the number of days between mixing events
$n_{pers}$ , the number of mixing periods to simulate
$n_{record}$ , the number of data points to record during each period. If the period plots look unsmooth, try increasing this parameter.

It is a good strategy to start with these default values, to see their results, then modify one of more to see the effects of those changes.

When you have the parameters you want, click the Run simulation button below to simulate the population and plot the results.

# Set up textboxes for simulation parameters
global t_mix_,n_pers_,n_record_
def solve_pars(t_mix,n_pers,n_record):
    global t_mix_,n_pers_,n_record_
    # Record ODE solving parameters for later use; this is executed each time
    # a parameter is changed
    t_mix_ = t_mix
    n_pers_ = n_pers
    n_record_ = n_record
    
out4 = widgets.interactive_output(solve_pars,{'t_mix':t_mix,'n_pers':n_pers,'n_record':n_record})
display(uis4,out4)

# Set up a button to run simulation with current parameters
buttonR = widgets.Button(
    value=False,
    description='Run simulation',
    disabled=False,
    button_style='', # 'success', 'info', 'warning', 'danger' or ''
    tooltip='Click to run the simulation from the current state, using the current parameters.',
    icon='' # (FontAwesome names without the `fa-` prefix)
)
outputR = widgets.Output()
@outputR.capture()
def solve(b):
    global D,t_mix_,n_pers_,n_record_
    # clear previous output
    outputR.clear_output()
    # Seed selected size class(es); isizes==[] means all size classes
    print('Running simulation...')
    D.solveODEs(t_mix=t_mix_,n_pers=n_pers_,n_record=n_record_)
    print('...done.')
    # Display output figures
    D.fig1.show()
    D.fig2.show()
    # Save images to files in png and svg formats
    D.fig1.savefig('period.png')
    D.fig1.savefig('period.svg')
    D.fig2.savefig('averaged.png')
    D.fig2.savefig('averaged.svg')
    #Image('period.png')

buttonR.on_click(solve)
widgets.VBox([buttonR,outputR])

References¶

Litchman, E., Klausmeier, C. A., & Yoshiyama, K. (2009). Contrasting size evolution in marine and freshwater diatoms. Proceedings of the National Academy of Sciences, 106(8), 2665–2670. 10.1073/pnas.0810891106